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Abstract. In this paper we examine how the predictions of conformal invariance can be widely exploited 
to overcome the difficulties of the density-matrix renormalization group near quantum critical points. 
The main idea is to match the set of low-lying energy levels of the lattice Hamiltonian, as a function of 
the system's size, with the spectrum expected for a given conformal field theory in two dimensions. As in 
previous studies this procedure requires an accurate targeting of various excited states. Here we discuss how 
this can be achieved within the DMRG algorithm by means of the recently proposed Thick-restart Lanczos 
method. As a nontrivial benchmark we use an anisotropic spin-1 Hamiltonian with special attention to the 
transitions from the Haldane phase. Nonetheless, we think that this procedure could be generally valid in 
the study of quantum critical phenomena. 
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1 Outline and General Facts 

The density-matrix renormalization group (DMRG) was 
invented by S. R. White in the early 90's and nowadays 
is recognized as one of the most accurate and efficient 
numerical techniques in the study of correlated quantum 
systems £Q. A number of authors is still contributing to 
its development for new applications and one of the ma- 
jor points is the dimensionality of the system. While for 
ID lattices of spins or electrons there exist several firm 
results, a generally accepted DMRG scheme in 2D is still 
lacking despite numerous proposals. However, even for ID 
systems, it is known that the DMRG encounters some 
difficulties in reproducing the physics of quantum criti- 
cal points. This is because one wants, at the same time, 
to explore larger and larger system sizes and to maintain 
a good numerical precision in order to locate the points 
where the excitation gap closes and to compute the as- 
sociated critical exponents. If the quantum chain has q 
states per site, the dimension of the full Hilbert space for 
L sites grows exponentially as q . The main approxima- 
tion of the DMRG is the truncation to the space spanned 
by the M eigenvectors of the block density matrix, p^, 
corresponding to the M largest eigenvalues. The problem 
is that, close to criticality, the system experiences fluctua- 
tions on very large scales and (as will be discussed in sec. 
0J) the spectrum of ph decays more slowly than in noncriti- 



cal cases. Hence if one considers, as a first indicator of the 
accuracy of the method, the sum of the discarded weights 
W M = £, > m w 3 ' then this different decay implies that 
the minimal M to reach a prescribed threshold is gener- 
ally larger when the system is close to a critical point. To 
be more specific, Legeza and co-workers found that 
the error in the ground state (GS) energy due to a finite 
M, at a fixed L, is simply: 

6E$ s (M)=E2 s (M)-E2 s (oo)<xW M . (1) 

Now, if we admit that the precision for the first excited 
state follows a similar trend, then there will be a regime, 
close to the critical point and/or for large L (see sec.0), 
where Wm becomes larger than the excitation gap itself. 
In addition, figs. 3 and 5 of ref. [2] indicate two other im- 
portant features (at least for the Ising model in a trans- 
verse field). First, at criticality the error in the energy 
at fixed M increases appreciably with the chain's length, 
and again we see that fixing M once for all, while taking 
larger and larger values of L, is not a completely safe pro- 
cedure. The second, and more important point is that the 
proportionality in eq. Q holds only after a few finite- 
system iterations have been performed. Using Legeza's 
terms, these are necessary to cancel the so-called envi- 
ronment error that dominates at sufficiently small L. The 
crossover from the environment-dominated regime to the 
truncation-dominated one is marked by a characteristic 
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size, L*(M), which increases with increasing M. In sec. El 
we will argue that a similar characteristic length emerges 
in critical spin-1 chains too. Moreover, in our opinion, the 
results of Andersson, Boman & Ostlund |2] regarding a 
critical chain of free spinless fermions, or equivalently a 
spin-1/2 XX model, point somehow in the same direction. 
Even if the system is known to be rigorously critical, the 
effect of a finite number of DMRG states is to introduce a 
fake correlation length, that grows as M 1 ' 3 . The interplay 
between an analogous DMRG length and the true corre- 
lation length near the critical temperature of classical 2D 
systems is also discussed in 0. 

We believe that the considerations above indicate that 
the critical behavior of (ID) quantum systems cannot be 
"simply" approached by using the computational resources 
to reach larger and larger values of L by means of the 
infinite-system DMRG, as it is sometimes done (see, for 
instance, and refs. therein). Rather, we prefer to ex- 
ploit as much as possible the consequences of conformal 
invariance. In many cases, indeed, we expect that the low- 
energy physics of a ID quantum system near its critical 
points is described by a suitable conformal field theory 
(CFT) in (1+1) dimensions |7ll%H51Tl7)l . Then, if the sys- 
tem under study has to be critical, we imagine that there 
will be a scale-invariant effective continuum model that 
captures its universal features. In 2D the key point is that 
all these universality classes are encompassed in the frame- 
work of CFT. In the last twenty years the numerous exact 
results in this area have been organized into an elegant 
framework, which is however too vast for our purposes. 
Therefore, in sec. 0we sketch only the points of contact 
with our analysis and refer to 1 1 1 j for a comprehensive 
guide and to [T^J for a shorter review on finite-size effects. 

However, for a correct interpretation of the DMRG 
data, and a comparison with a suitable CFT, an accu- 
rate calculation of the first excited states at varying L is 
needed. The targeting of more than one state was pro- 
posed by White himself JHj and used by different authors 
in subsequent papers. A recent comment on this prob- 
lem can be found in ^3] • O ne of the aims of our work is to 
test, in connection with the DMRG, the recently proposed 
Thick-restart Lanczos method £01 to handle a sufficient 
number of excited states. 



2 Identification of the CFT through the 
Finite-Size Spectrum 

The requirement of scale invariance in 2D is sufficiently 
strong to allow (under general assumptions) a classifica- 
tion of the states of a quantum field theory in terms of the 
irreducible representations of the Virasoro algebra gener- 
ated by a set of operators satisfying: 



[L m ,L n ] = (m-n)L m+n + —5 m+nS ){m z -m) , m,nEZ 

• (2) 

(As is customary in CFT, the statements and equations 
regarding the holomorphic part should always be suitably 



repeated for the antiholomorphic part, denoted by over- 
bars). For unitary theories the central charge of the alge- 
bra is either c > 1, or can be chosen as one of the values 
c = 1 — 6/p(p +1) (p = 2, 3, . . . ), each one corresponding 
to a so-called minimal model. In the latter case, once a 
c < 1 is given, one can decompose the Hilbert space into a 
finite number of irreducible representations of (J2J labeled 
by the eigenvalues of the generator Lq: 



Lq\A) = A\A) , A 



[jp + l)r-ps] 2 - 1 
MP + 1) 



(3) 



with l<s<r<p — 1 € Z. More specifically, the pri- 
mary states |Z\) are annihilated by L m with positive m 
and each one of them generates a conformal family that 
contains all the states obtained through the application of 
the negative-integer generators \A) k n = (L^ m ) k \A) (the 
secondary or descendant states). It can be seen, as a con- 
sequence of the Virasoro algebra, that these are in turn 
eigenvectors of Lq: 



L \A) k m = (A + mk)\A) k m 



(4) 



In principle, all the correlation functions of the theory can 
be reduced to correlators of primary operators. Moreover, 
the operators Lq and Lq play a special physical role as 
can be seen by considering a conformal transformation 
that maps the plane onto an infinite cylinder, whose cir- 
cumference of length L represents the space axis with pe- 
riodic boundary conditions (PBC). Then the energy and 
momentum operators are simply expressed as: 



Hcft = 



2ttv ' 
~L~ 



Lq 



12 



2tt 
~L 



[Lq 



Lq] . 

(5) 

Actually, the pre-factor v has been inserted "by hand" 
in view of the identification of the (critical part of the) 
spectrum of the original quantum Hamiltonian with the 
values predicted by the CFT construction. More specifi- 
cally, the ground state is simply identified with |vac), for 
which both Lq and Lq have null eigenvalues. Thus, the 
finite-size corrections to the vacuum energy are: 



7TCV 

~6L 



(6) 



whereas the excited states follow from the construction of 
eqs. © and |TJJ): 

E L (A,m,k;A,fh,k)-E2 S = ^-[A+A+mk+fnk] . (7) 

L 

It is important to recall that the sum of the eigenvalues 
of Lq and Lq, that is the last term in square brackets 
do = (A + A + mk + jfifc), is called the scaling dimension 
because it enters the algebraic decay of the correlation 
function of the corresponding operator O: 



-2do 



z—z—r 



(O(0,0)O(z,z)) ~ z -nA+mk)z-2{A+mk)\ 

(8) 

In the language of renormalization group (RG) theory, the 
operators are seen to be relevant when do < 2, irrelevant 
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when do > 2 and marginal when the scaling dimension is 
exactly 2 (in 2D). The exponents, y*, of the scaling fields 
near a fixed point of the RG flow are simply given by 
y. L = 2- di. 

The edge case c = 1 opens the way to unitary CFT 
with an infinite number of primary operators. It is re- 
lated to certain topological constructions of a free bosonic 
field ^J, and the operators (or the associated states) can 
be further classified in terms of an U(l) Kac-Moody al- 
gebra. Here we won't enter into more details since ref. 
[*T7) is specifically devoted to the study of such c = I 
CFT starting from the spin-I Hamiltonian i|12|) . using just 
the method explained in the present paper. Finally, when 
c > I we have a continuum of allowed values and the com- 
plexity of the problem in its generality simply forbids us 
to go on with the discussion of classified cases. We refer 
to the literature on CFT for general discussions (see, for 
example, (lllllfcip . We will content ourselves with the ob- 
servation that, at least from an operative point of view, in 
many cases eqs. ©-(jSJl should remain valid once the labels 
A, A and the secondary indices are properly interpreted 
in terms of wider symmetries that have to be investigated 
case by case. 

The first problem is, of course, that one does not know, 
a priori, which is the CFT that has to be invoked. A great 
step towards the answer to this question is made if one 
knows the central charge of the theory, c, that character- 
izes the underlying Virasoro algebra. Apart from eq. J2J) 
and the operator product expansion of the stress-energy 
tensor, the central charge appears explicitly in the size- 
dependence of the GS energy density: 



E2 S 



TTCV 

6L2 



(9) 



which is nothing but eq. © plus an infinite-size term, 
which is absent in the formal CFT but has a definite 
value for a given quantum Hamiltonian. If the quantum 
chain has open boundary conditions, which are generally 
better for DMRG convergence, eq. © should be modi- 
fied in two ways 18,19 . The denominator becomes 24L 2 
and a boundary term B/L should be added, with a non- 
universal pre-factor B. This term introduces a slower con- 
vergence of the GS to the thermodynamic limit (TL) and, 
at least in this sense, PBC may be more useful for finite- 
size scaling (FSS) [201 - 

In quantum field theory the second term of @ ac- 
counts for the Casimir effect [T^l, while in statistical sys- 
tems it gives the correction to the free energy at small 
temperatures (L playing the role of 1/T [2J)- In our case, 
eq. © is the starting point to discover the CFT that is 
appropriate for the problem under consideration. Indeed, 
having good estimates of _E^ S at various L, we can first 
best-fit eoo and cv. Then, if we have reasons to believe 
that one of the states has scaling dimension 1 (eqs. J7J 
and ©), then we can calculate: 



If we now interpret Ak = 2tt/L as the quantum of mo- 
mentum for a chain of length L, then eq. (|10|) looks like 
a discrete derivative of the energy vs momentum relation, 
that is, a group velocity. This somehow explains the term 
"spin velocity" , which is widely used in the literature of 
quantum spin chains independently of the real nature of 
the excitations. Again, the exact form of eq. 1)10(1 may vary 
depending on the boundary conditions. In spin chains with 
twisted boundary conditions the general relation in eq. (J7J) 
remains valid and it's the scaling dimensions of the oper- 
ators that turn out to depend on the twisting angle \22\. 
With open boundary conditions, according to ref. JB] the 
essential modification is the replacement of L by 2L. This 
justifies the coefficient of L 2 in eq. |j5Jl and eq. (58) of |19j . 
that looks like eq. 1(10(1 with a quantum of momentum 
Ak = tt/L. 

Summing up, with a combined usage of eqs. (JHJ and 
(110(1 we can calculate c from the numerical data of the 
first excited levels for different chain's lengths. Then the 
identification may fall into two cases. For c < 1, unitarity 
demands that the values of the central charge are quan- 
tized, according to the list of minimal models. In this case, 
a small discrepancy from one of these values is likely to 
be related to numerical uncertainties. The matter becomes 
more complicated when one finds a numerical c larger than 
one, in which case unitarity by itself is not enough to pro- 
vide a quantization condition. The symmetries of the lat- 
tice Hamiltonian in this case are of great help because we 
expect them to be present also in the corresponding con- 
tinuum model. Then one may focus on the known 2D field 
theories whose actions are invariant under both this sym- 
metry group and under conformal transformations. Be- 
side that, when one studies a novel and computationally 
demanding system it may be very useful to get a first in- 
sight of the CFT by performing a finite-size analysis of 
DMRG data obtained with the infinite-system algorithm 
and open boundary conditions, as in [T2|. Nonetheless, we 
discourage a blind usage of the infinite-system DMRG to 
conclude (as in ref. [H]) that the numerical discrepancies 
from the expected values are to be ascribed to nontrivial 
effects beyond the CFT framework. If in doubt, the final 
answer should always come from a careful refinement us- 
ing finite-system numerical data, in a range of L where 
the scaling behavior is visible but the accuracy is not cor- 
rupted significantly by the DMRG truncation error. 

In any case, what we are actually performing is a self- 
consistent guess of the underlying CFT. In principle, this 
allows an analytical calculation of the scaling dimensions, 
di : of all the operators associated with the excited states. 
The corresponding energy gaps, here formally indexed by 
£, scale according to eq. (JJJ): 



AE\ 



2tt 



-vd e 



(11) 



AE L {d=\) 
2ir/L 



(10) 



Matching the numerical spectrum of a certain number of 
levels with the structure encoded in the scaling dimensions 
of eq. ((11(1 represents a particularly stringent test of the 
hypothesis made. 
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3 Model and Details of the Implementation: 
Multi- target Method 

Following our previous paper |17j we have considered the 
following spin-1 Hamiltonian: 



H 



E 

3 



(sts- 



3+1 



Sj s j+lj 



D{S 



z\2 



(12) 

for a chain of L spins which includes both an Ising-like 
and a single-ion anisotropy term, with coefficients A and 
D respectively. Formally we have set h — 1 and the overall 
coupling constant J = 1 so that every quantity turns out 
to be dimensionless and we have imposed PBC: Sj+l — 

for a discussion 



of the various phases in the X-D diagram. 

The algorithm that we have implemented for DMRG 
calculations follows rather closely the lines indicated by 
White in his seminal papers [T31I25] . As regards our spe- 
cific application, we should outline the following points. 

The superblock geometry was chosen to be [B r • | B^ ef •] 
with PBC, where B^ of is the (left <-> right) reflection of 
block B r with r' sites. The rationale for adopting this 
configuration is that, being effectively on a ring, the two 
blocks are always separated by a single site, for which 
the operators are small matrices that are treated exactly 
(no truncation) I13| . In this way we expect a better preci- 
sion in the correlation functions calculated by fixing one 
of the two points on these sites and moving the other one 
along the block. In a recent application of the DMRG to 
quantum chemistry calculations [Sj it has been pointed 
out that the configuration of the superblock may be one 
of the major points of optimization of the method for fu- 
ture applications. As regards the choice of the boundary 
conditions, we are aware of the fact that with open con- 
ditions a smaller M is generally required and that in cer- 
tain cases (i.e. Gaussian transitions) the introduction of 
twisted boundary conditions is a clever trick to identify 
the critical point using numerical data on relatively small 
systems (221127ll2l31l37)] . Nevertheless, in our set of calcula- 
tions we have adopted PBC to get rid of the edge effects 
that in some cases mask almost all the information con- 
tained in very short-ranged (string) correlation functions. 
Moreover, we are primarily interested in the transitions 
from the Haldane phase, for which it is known that the 
finite-size GS acquires a fourfold degeneracy due to the 
two free effective spins at the ends of the chain . With 
PBC instead, one has a unique finite-size GS so that the 
convergence of the Lanczos procedure is better and the 
analysis of the numerical gaps is simpler. 

We used the finite-system algorithm with three iter- 
ations. This prescription should ensure the virtual elim- 
ination of the so-called environment error 0, which is 
expected to dominate in the very first iterations for L < 
L*(M) (see below). In fact, it is only after a suitable num- 
ber of such sweeps that we may expect that the error in the 
energies has been minimized (for a given L and M) and eq. 

holds true. Normally the correlations are computed at 



the end of the third iteration, once the best approximation 
of the GS is available. 

Dealing with quantum spin chains, we always exploit 
the conservation of the ^-component of the total spin, 
M z . Typically we are interested in nonmagnetic GS's, 
that is, with M z = 0. In every studied case the corre- 
lations have been calculated targeting only the lowest- 
energy state within this sector. However, in order to an- 
alyze the energy spectrum, we had to target the lowest- 
energy state(s) in the other sectors \M Z \ = 1,2,... and/or 
target also a few excited states within the M z = sec- 
tor, depending on the phase under study. The standard 
Lanczos algorithm gives with enough precision the ground 
state of the system, but it is not so accurate for the excited 
states. Moreover, as a consequence of a general theorem on 
tridiagonal symmetric matrices, we can't have degenerate 
states from this method. So the algorithm must be mod- 
ified to allow the building of the reduced density matrix 
over the block as a mixture of the matrices correspond- 
ing to each target state. While for the latter point we 
are not aware of any specific "recipe" other than that of 
equal weights, the implementation of the multi-target di- 
agonalization routine within our DMRG code is based on 
the Thick-restart Lanczos method recently introduced by 
Wu and Simon |15| . In principle, the strategy used here 
does not rely on a particular algorithm to extract a group 
of eigenvalues of the superblock Hamiltonian. The meth- 
ods commonly used so far are the Davidson-Liu and the 
block Lanczos. Our choice of the Thick-restart Lanczos 
was based on the intention of testing and exploiting the 
stability of the method, as claimed by the authors |15j . 
Moreover, its implementation is a simple extension of the 
standard Lanczos procedure, as described in appendix A. 

Once M z is fixed, in a given run we wish to follow 
simultaneously the first levels \M z ;b) with b=0,l,2,. . . ,t 
(the GS being identified by (M z = 0;b= 0)). Then, as 
in the conventional Lanczos scheme, we have to iterate 
until the norms of the residual vectors and/or the differ- 
ences of the energies in consecutive steps are smaller than 
prescribed tolerances (10 -9 — 10 -12 in our calculations). 
The delicate point to keep under control is that, once the 
lowest state \M Z ;0) is found, if we keep iterating search- 
ing for higher levels the orthogonality of the basis may be 
lost, just because the eigenvectors corresponding to these 
levels tend to overlap again with the vector |M Z ; 0). As a 
result, the procedure is computationally more demanding 
because one has to re-orthogonalize the basis from time to 
time. We have seen that this part takes a 10-20% of the 
total time spent in each call to the Lanczos routine. We 
have also observed that if this re-orthogonalization is not 
performed, one of the undesired effects is that the excited 
doublets (typically due to momentum degeneracy) are not 
computed correctly. More specifically, it seems that while 
the two energy values are nearly the same in the asym- 
metric stages of the sweeps, when the superblock geom- 
etry becomes symmetric (r = r' in the notations of the 
preceding point) the double degeneracy is suddenly lost 
in a spurious way. This is in line with the results of ref. 
Pj, where the error in the energy is kept under control 
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by means of a dynamically-adjusted M and the configu- 
rations that require a larger number of DMRG states are 
just those near the symmetric one. 



4 Dependence on the Number of DMRG 
States 

As discussed in the introduction of sec. 0] the choice of 
the number of optimized states with respect to the chain 
length L is the crucial point to address in any DMRG 
calculation. Being conscious that the energy accuracy is 
not an exhaustive indicator, from eq. Q onc should try, 
at least, to keep the discarded weight Wm as small as 
possible. This quantity, in turn, is related to the decay of 
the density matrix eigenvalues {uij} as a function of the 
index j: 

j>M j 

where the last equality is a simple rewriting in terms of 
the degeneracies Dj and of the distinct eigenvalues Zj . The 
issue appeared to be important from the very beginning 
and is the core to the success of the DMRG. In |13U31| 
it was argued that the convergence of the GS energy of a 
gapped or a spatially finite system is roughly exponential 
in M : 

SE L oc e - M / M *W , (14) 

with a superimposed step-like behavior, probably related 
to the successive inclusion of more and more complete 
spin sectors. It is also generally believed that this expo- 
nential decay becomes slower (possibly algebraic) when a 
critical point is approached. At this stage it is interest- 
ing to recall that for integrable systems the spectrum of 
Pb can be determined exactly. In fact, it is known |5f2| 
that for a quantum chain the infinite-size block density 
matrix is given by pb — X 4 , where \ is the corner transfer 
matrix of the associated 2D classical statistical system. 
This statement has a wide generality, with the exception 
of critical cases where the boundary effects may have some 
role and are expected to affect the tails of the distribution 
|32|. For integrable systems a further step can be made: x 
is expressed as the exponential of a pseudo-Hamiltonian, 
K, that involves the same local operators of the Hamilto- 
nian (e.g. S"S" +1 ) but with coefficients depending on the 
site index j. Moreover, the spectrum of K can be deter- 
mined exactly and, typically, the eigenvalues turn out to 
be equally spaced. Therefore the distinct eigenvalues of pb 
decay as: 

Zj oc Z 3 , Z = e" e , (15) 

e/4 being the level spacing of — A'. These predictions have 
been explicitly verified for the Ising model in a transverse 
field (Hj and for the XXZ spin- 1/2 Heisenberg chain PUT] . 

Now we shall try to elucidate this topic in the case of 
critical spin-1 systems using the Hamiltonian of eq. (|12(l . 
At (A = 1, D — 0) one has an isotropic anti- ferromagnetic 
(AF) Heisenberg chain of integer spin, whose theoreti- 
cal interest comes from what is now known as Haldane's 



conjecture |33| . that predicts a genuine quantum behav- 
ior with a finite energy gap in the excitation spectrum. 
This has to be contrasted with what happens in half-odd- 
integer cases that are gapless (i.e. critical) in the ther- 
modynamic limit. Another important theoretical contri- 
bution is the mapping between spin chains and restricted 
solid-on-solid models proposed by den Nijs and Rommelse 
|34| . Qualitatively, the Haldane phase is interpreted as a 
spin-liquid, in the sense that the effective particles - the 
spin state |0) represent an empty site and |±) represent 
a particle with spin up or down - are positionally disor- 
dered but carry anti-ferromagnetic order. In the quantum 
states that contribute to the GS, this order is hidden by 
arbitrarily long strings of |0)'s but can be measured by 
the string order correlators j26U34| : 

0§(j, k) = - (sf exp ^itt £ S°\ S%j , (16) 

(the expectation value being taken on the GS) and their 
asymptotic order parameters, Oj? = lim|j_fc|_ >0O Og (j, k). 

Beside that, nowadays there exist several experimen- 
tal realizations of these anisotropic spin-1 chains: to the 
author's knowledge, RbNiCl 3 |2BJ and Y 2 BaM0 5 [351136] 
are examples of pure Haldane systems, CsNiFe3 is a ferro- 
magnet (A = —1) with appreciable single-ion anisotropy 
(D ~ 0.4) and C0CI2 2 H2O behaves as an Ising ferromag- 
net (|A| ^> 1) with high easy-axis anisotropy (D ~ —5) 
[33137] . The so-called NENP [331135] and NENC [UJ rep- 
resent, respectively, small-D (~ 0.2) and large- 1? (~ 7.5) 
antiferromagnets with easy-plane anisotropy. 

In order to study the convergence of truncation errors 
close to criticality we have selected a pair of representative 
points in the phase diagram. In the origin (A = 0, D — 0), 
where a Berezinskii-Kosterlitz-Thouless (BKT) transition 
is expected |3T], we have the spin-1 analogue of the XX 
spin- 1/2 model used in ref. fj] to examine the effect of a 
finite M on a quantum critical system. The second point 
that we will examine in this section is (A = 1,D = 0.95) 
because from refs. [171130] we know that it lies very close 
to the line of transition from the Haldane phase to the 
so-called large- D phase. In the following numerical anal- 
ysis these two points will be denoted by XX and H-D, 
respectively. The results can be thought as the critical 
counterpart of the ones presented in |10] for an isotropic 
spin-1 Heisenberg chain. 

For the H-D point we have computed the first excited 
state in the sectors with M z = 0, 1 for L = 16, 20, 24, 32, 48 
and 64, using M = 81, 162, 243, 324 and 405 DMRG states 
for each case. For the XX point, we have monitored the 
same states for the same values of M using also L = 28. 
The exponential decay of eq. ifHjl seems to be adequate 
in all the cases, as shown in fig. ^ For every fixed L we 
best-fit the energy-vs-M data and read off the character- 
istic values M*{L), represented in fig. [21 Of course, the 
larger L is the larger M* is, and it seems that the curves 
do not saturate indicating, as is reasonable, that a finite 
M is not enough to describe properly a (quasi-)critical 
system with arbitrarily large L. However, the promising 
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Table 1. Slopes Mo of the best fits of the sets plotted in fig. 
[5] (with same notations). ' Only with L > 28; J Only with 
L > 24. 



Data Set 




Mo 


GS XX 




16.9 ± 0.4 


1st M z = 





17.6 ± 0.4 


1st M z = 


1 t 


15.6 ± 0.4 


GS H-D 




13.2 ± 0.3 


1st M z = 





15.0 ± 0.2 


1st M z = 


i i 


11.4 ± 0.3 




M 



Fig. 1. Energies of the GS and of the first excited states within 
M z = 0, 1 at the points XX and H-D (see text). The symbols 
represent DMRG values obtained with L — 32 and an increas- 
ing number of DMRG states while the continuous lines are 
exponential fits. For clarity reasons, the offsets reported in the 
legend have been added. 
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Convergence NOT reached 








i 



ln(L) 



Fig. 2. Exponential factor M*(L) (eq. 1141 1 for the points 
indicated in the legend (see text). The lines separate two re- 
gions where the convergence of the low-lying levels with M is 
essentially reached or not. 



feature of fig.[21is that the various data sets approximately 
lie in a strip of the M — lnL plot, characterized by a lin- 
ear slope Mq = 15 ± 3. More precisely, this value is an 
overall measure of the slopes of the best-fit straight lines, 
reported in table for the six sets considered. In other 
terms, if we invert the relation between M and L, we find 
L*(M) oc exp (M/M ), which means that for fixed M we 
expect a good convergence of the energies for L < L*(M) 
(negligible truncation error, in the language of 0) and 
that the gain in increasing M is exponential with a sur- 
prisingly small reference value Mo (in line with the first 
observations of White himself Interestingly enough, 
the logarithmic behavior of M* (L) can be justified by tak- 
ing eq. (|15J) to be valid even at critical points. There one 
expects e — > and indeed conformal invariance indicates 
that the finite-size level spacing vanishes as e = eq / In L 
(with eo a constant) jjjj. Now, if we approximate the dis- 
carded weight of eq. (|T3)l ignoring the degeneracy factor 
Dj, we readily get: 

W M = Y = —Z M ~ ^L e -^M/lnL (17) 

1 — Z en 

3>M 

Unfortunately, the effective accuracy gets poorer, by 
one or two orders of magnitude [3J, when we deal with 
correlation functions. Keeping the errors in the low-lying 
levels below the desired threshold may not be sufficient 
so that we had to adopt an additional criterion to decide 
whether the selected M is large enough or not. In practice, 
we check systematically the properties of translational and 
reflectional invariance that we expect from the symme- 
tries of the Hamiltonian. In fact we have observed that 
of one the "symptoms" for a too-small M is the visible 
(i.e. above numerical uncertainties) lack of some of these 
invariances. To be more specific, if C(0, k) is a certain 
correlation function computed starting at j = 0, we have 
always increased M (at the expense of L) until the bound 
\C(L/2,L/2±k) - C(0,fc)|/|C(0,fc)| < 0.05 was met for 
k varying from to L/2, possibly with the exception of 
the ranges where C(0, k) is very small, say 10 -6 . In fig. El 
we give an example with string correlation functions (eq. 
(I16fl 1 at the XX point. These should be translationally in- 
variant but their numerical estimates depend in fact on 
the starting point j because we have intentionally fixed a 
too-small value, M = 50. It is also interesting to point out 
that in this example Og(j, j + r) essentially coincides with 
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(— ) r (SjSj +r ) (not plotted). We interpret this coincidence 
with the onset of planar order at the BKT transition. At 
the XX point both the transverse and the longitudinal 
string order parameters are expected to vanish (sec. II of 
|42p. and the fact that the minimum of the x-data in fig. 
131 is about ~ 0.3 is to be interpreted as a finite-size effect. 
In a similar way, a preliminary study of the ID Hubbard 
model with bond-charge interaction near its critical points 
03] indicates that in order to obtain accurate correlation 
functions it is necessary to use a large number of DMRG 
states (> 500) and several finite-system iterations. From 
the literature it turns out that this is a general feature of 
electronic systems, like Hubbard models and critical vari- 
ants. 

Finally, we mention the relevance for ref. 04 , where 
the thermal behavior of the quantum ID and 2D S = 1 
XX model is studied by means of the two-time Green's 
function method. The adopted decoupling scheme requires 
the calculation of the on-site, nearest-neighbor and next- 
to-nearest-neighbor ordinary correlations in the x and z 
channel: 



o = 



(18) 

(as in eqs. (14) and (15) of 03] even if it is not clear if Cf 
should include the on-site term with S = —5' or not). It is 
explicitly said that an exact term of comparison in ID at 
T = would be needed in order to check more precisely 
their methodology. Though in principle approximate, the 
DMRG data can be regarded as a solid benchmark. We 
have computed the five numbers above with M — 300 
for L = 32, 48, 64, 80, 100. Having to do with short-range 
non- vanishing quantities, the size dependence is very weak 
(fourth decimal place or better) and an algebraic extrap- 
olation to L -► oo yields C$ = 0.7577, C$ = 0.4846, 
Cf = 0.5579, CI = -0.1778 and Cf = 0.464 (not in- 
cluding the on-site term). From a direct comparison with 
table I of 03] we see that with their choices of the decou- 
pling parameters the largest discrepancy affects Cf with 
a relative difference of about 20%. 



5 Examples and Results 

The quality of the numerical analysis of the critical prop- 
erties depends heavily on the location of the critical points 
of interest. At present, our study focuses primarily on the 
transitions from the Haldane phase, for which it is con- 
venient to fix some representative values of A and let D 
vary across the phase boundaries. This preliminary task 
of finding D C (X) turns out to be crucial for subsequent 
calculations and is divided into two steps. 

First, one gets an approximate idea of the transition 
points using a direct extrapolation in l/L of the numer- 
ical values of the gaps, computed at increasing L with a 
moderate number of DMRG states. Clearly, one may want 
to explore a rather large interval of values and so the in- 
crements in D will not be particularly small (say 0.1). An 




Fig. 3. Transverse (a = x, circles) and longitudinal (a = z, 
squares) string correlation functions for a chain of 64 sites with 
PBC at the XX point (see text), computed with only M = 
50 DMRG states. The empty symbols represent Og(j = 0, r) 
while the full ones represent Og (j — 31, 31+r), with an evident 
dependence on j. 



0.25 - 



— AT 


= 


-- M~ 


= 1 


■■■■ M z 


= 2 


— M z 


= 3 



3 



0.15 




Fig. 4. Extrapolations of the gaps between the GS and the first 
excited states within the indicated sectors of M z at A = 0.5. 
The values for L — > oo have been obtained through a quadratic 
best fit in l/L from the data at L = 32,48,64,80 with 216 
DMRG states. For graphical convenience, the various sets have 
been turned into continuous curves using splines. The vertical 
gray line indicates the point of intersection between the states 
with M z — and M z = 1, that is, the analogue of the Haldane 
triplet at A = 0.5. 



example of such a scanning at A = 0.5 is presented in fig. 
0] Note that at (A = 0.5, D ~ 0.6) one has a first insight 
of the "cascade" of levels predicted by CFT (eq. Q for 
L — » oo.) 

As a second step, the analysis must be refined around 
the minima of the curves AE-vs-D with smaller incre- 
ments in D and a larger value of M. According to FSS 
theory ( |45| and appendix B), at the true critical point 
(that is, in the TL) one should see that the data settle 
to a constant in the log-log plot of the scaled gaps vs L. 
Fig. [S] shows an example of such an inspection for A = 0.5 
and D varying about the H-D transition. It is seen that 
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ML) 

Fig. 5. H-D transition at A = 0.5: Scaled gaps (between 
\M Z = 0,b = 0) and \M Z = l,b = 0» at L = 
10, 12, 14, 16, 18, 20, 22, 24, 32, 48, 50 (with 400 DMRG states) , 
for the values of D indicated in the legend. 



the differences in the slopes of the various curves are not 
so pronounced. Hence, from this plot we have selected two 
candidates for the critical point _D C (0.5), namely D = 0.62 
and D = 0.65. At this stage we should mention that in 
this type of transition the phenomenological renormal- 
ization group (PRG) method 45J also typically yields a 
pair of (pseudo-) critical values at fixed L. In this case, 
these two sequences of values seem to converge to the 
point D = 0.62 (but in order to see convergence of the 
curves up to L = 50 we had to use 400 DMRG states). 
Unfortunately, this turns out to be an invalid tie-break 
because at this point the finite-size (3- function (eq. J33)) 
increases with increasing L. This is shown in fig. El where 
/3l(0.65) is also plotted. The latter scales to zero with a 
size-dependent slope that we calculate from eq. I|38|) . The 
extrapolation to 1/L — > (and restricted to L > 22) then 
gives v = 3.69 ± 0.04. As discussed below this is not a 
particularly good estimate of v. Nonetheless, the location 
of the critical point D c (0.5) = 0.65 is quite close to the 
value D = 0.635 obtained in ref. [21] with a method based 
on twisted boundary conditions and exact diagonalization 
up to L = 16. In this sense, we suspect that the difficulties 
encountered both with FSS and PRG (roughly speaking, 
the "splitting" of critical points) are due to the peculiar 
structure of the energy spectrum in this type of transition. 
In particular, both sides of the transition are massive and 
it is likely that with PBC we are faced with the scenario 
proposed by Kitazawa [22]: In eq. (j3l)]) the constant 
for the first excited state could vanish and we have to 
consider a second-order expansion in (D — D c ). Conse- 
quently, the values of the critical points are determined 
via parabolic intersections that are more sensitive to nu- 
merical uncertainties. Probably this is also the reason why 
the scaling analysis of (3l as used here performs poorly. 
In fact, in our framework the best estimate of v comes 
from another method, namely via the scaling dimension 
of the mass-generating operator. According to the discus- 




ln(L) 

Fig. 6. Finite-size /3-functions for (A = 0.5, D = 0.62) (open 
squares), (A = 0.5, D = 0.65) (full squares), (A = 0.5, D — 
— 1.2) (open circles), (A = 1,D — —0.315) (full circles) and 
(A = 1.185, D = 0) (full diamonds). All the numerical deriva- 
tives have been calculated by means of centered differences 
with 5D = 0.01 for the first two cases (Haldane-large-D tran- 
sition), SD = 0.005 for the second two and S\ — 0.005 for the 
last one. Notice that the three lines for the transitions from the 
Haldane to the Neel-like phase (Ising type) have essentially the 
same slope, that is, the same v (see also the discussion at the 
end of sec. [SJ. 



sion that follows, for (A = 0.5, D — 0.65) we find v = 2.38, 
essentially in accord with the values given in |24| . 

At the transition between the Haldane and the Neel- 
like phase (henceforth denoted as H-N) these difficulties 
regarding the location of the critical points and the scaling 
to zero of (3l are not experienced. Let us work out in detail 
the A-driven transition at D = 0, that has been recently 
revisited [S] to argue that it does not belong to the Ising 
universality class, as is generally accepted. Taking advan- 
tage of previous results, we fixed D = in eq. I|12|l and 
varied A about 1.18. In fig.Efa) we report the L-dependent 
pseudo-critical values obtained by solving numerically the 
PRG equation 45 : 



\(L + 6L)AE L+SL {\ pc ) - LAE L {X pc )} 
SL 



= 



(19) 



Our best-fit critical value in the TL turns out to be A C (D = 
0) = 1.1856 (see caption), in agreement both with |2"9~] 
and with 6 . This point separates a gapfull phase, where 
the scaled gap increases with increasing L, from a gapless 
one (doubly degenerate GS) , where it is expected that the 
scaled gap converges rapidly to zero [23]. The algebraic 
decay of /3l is rather evident from fig. Eland one can read- 
ily estimate v = 0.987 ± 0.002 from the slope of the linear 
best fit. This is the first indication that the H-N transition 
belongs to the (2D) Ising universality class, as reported 
by various authors (see below our argument based on the 
combined use of CFT and DMRG). 

Despite some difficulties in locating the critical points 
of the H-D line on the basis of DMRG data, we are now 
in position to suggest an extension in the usage of this mi- 
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Inverse size 

Fig. 7. (a): PRG location of the Ising-like critical point A c 
at D = 0. The pseudo-critical values of A (full circles) are 
plotted against the inverse of the middle size, (L + <5L/2) _1 , 
and correspond to L ranging from 12 to 48 in steps of SL — 4 
(see eq. 119H ) using M = 405. The continuous line is a best- 
fit of the form A c — ALcxp (— L/Z/prg) with A — 0.1185 and 
LpnG ~ 3.224. (b): Extrapolation (continuous line) of the GS 
energy density (full circles) according to the form discussed in 
the text. 



merical method to approach the critical points of quantum 
Hamiltonians. Taking advantage of the underlying con- 
formal theory that should be present in all these cases, 
we focus on the finite-size energy spectrum, as summa- 
rized by eqs. (0 and Once the critical point is lo- 
cated, we select a number of states that seem to become 
degenerate with the GS in the limit L — > oo. Then we 
consider LAEl/2it and plot the data against 1/L to see 
whether they settle to a constant, which represents the 
associated scaling dimension de multiplied by the veloc- 
ity v. Actually, due to this pre- factor we have to imagine 
a self-consistent procedure: Depending on the type of the 
transition we have in mind (that is, depending on the con- 



formal anomaly c), we stick on one or more levels in the 
spectrum that have exactly d = 1. Hence the value for 
1/L — * is nothing but v. Once the velocity is estimated, 
one uses eq. (0 to best fit the product cv and see if the 
value of c and the hypothesis concerning the universality 
class are self-consistent or not. 

To clarify the matter, let us return to the H-N transi- 
tion, that is thought to be in the universality class of the 
classical 2D Ising model with c = 1/2. This can be con- 
sidered as the paradigm of minimal models in CFT and 
the corresponding quantum field theory is that of a mass- 
less Majorana fcrmion. Apart from the identity (with zero 
scaling dimensions) one has only two primary operators 
with Ah = 1/16 and At — 1/2. In addition, modular in- 
variance 0S| demands that the conformal spin (A — A) 
of the combinations that enter the partition function (on 
the torus) must be an integer (zero in this case). Hence 
we are left with two nontrivial primary operators of scal- 
ing dimensions dr = 1 and dh = 1/8, that are inter- 
preted as energy and spin density, respectively. The for- 
mer is responsible for the variations away from the critical 
temperature and so the correlation length index is given 
by v = 1/(2 - dr) = 1. The latter, from eq. ©, gives 
the decay exponent of the spin-spin correlation function 
rj z = 1/4 at the critical point. At the same time, being 
v = 1 (see the following eqs. I|39|) and l|40[l in appendix 
B), one determines also the exponent j3 = dh = 1/8, de- 
scribing the opening of the magnetization away from crit- 
icality. 

The GS energies at various L for the critical point 
(A = 1.1856, D = 0) found above are plotted in fig. E^b). 
We have seen that in the range L = 12—56 the small devia- 
tions from —1.5 are accurately fitted by the form E2 S / L = 
e 00 -CT7r/6L 2 -A/ J L 3 / 2 exp(-L/^i), with e x = -1.50093, 
cv = 1.327, £i = 3.48 and A — 1.251. This dependence can 
be justified by recalling that the lattice model is mapped 
onto a CFT (yielding eq. ©) for the critical sector, plus a 
residual Hamiltonian that accounts for the massive levels. 
The exponential term is what one expects in a massive 
regime 07|, and we observe that the length £i multiplied 
by the finite gap with the sector \M Z = 1|, AE?W = 0.754, 
yields v\ — 2.62, a typical value of velocity (close to the 
one found below). In fact, at this stage we observe a non- 
trivial feature: the massless modes described by the CFT 
seem to be associated only with the levels within M z = 0, 
while those with M z ^ maintain a finite energy gap 
in the TL. Hence, the reference state for the calculation 
of v will be the second excited state in M z = 0, cor- 
responding to conformal dimensions (1/2,1/2). Using eq. 
(|10fl we determine a (CFT) velocity that extrapolates to 
v = 2.676 ± 0.001 (see fig. |EJl, and the resulting central 
charge, c = 0.4959 ± 0.0006, lead us again to the 2D Ising 
universality class as in 121] and refs. therein. An even more 
stringent confirmation comes from fig. |H| that illustrates 
the core of the multi-target method proposed here. In ta- 
ble [21 we compare the theoretical values with the numer- 
ical estimates of the scaling dimensions (second column) 
obtained by taking the ratios between the intercepts at 
1/L = of the data set in fig.|Hland the velocity pre- factor 
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v = 2.676. These intercepts, in turn, are determined using 
power-law fits to eliminate the residual dependence on L 
in a fashion similar to that of |30|. The fits are carried out 
in the range L > 36 because of the shoulder at L ~ 40 
in fig. [7| that signals the onset of the scaling regime. All 
the multiplicities are met, even if with DMRG calculations 
alone we are not able to classify the four degenerate states 
in terms of the secondary indices. The fourth column con- 
tains the total momentum/conformal spin expected from 
eq. Ipfl. Question marks indicate that the conformal con- 
tinuum theory predicts also for those cases that are ex- 
pected to have \Q\ — n We suspect that this is due 
to the correspondence between the original spin model 
and the field theory that maps the discrete AF structure 
(|Q| = tt) onto the low-momentum sector. Indeed, for the 
CFT to give \Q\ — tt a secondary index equal to L/2 
would be needed, and this would yield a nonzero energy 
gap in the TL. The overall agreement is good (1% in the 
worst cases) and the complete structure of the spectrum 
of the c = 1/2 minimal model is reproduced (only the rel- 
evant and marginal cases with d < 2 are reported). Note 
that all the marginal operators have nonzero momentum 
and so they cannot represent a valid perturbation to the 
continuum Hamiltonian because they would break transla- 
tional invariance. The absence of marginal operators sug- 
gests that each point of the H-N transition corresponds to 
the same c = 1/2 theory and the line in the phase dia- 
gram is "generated" by the mapping of the discrete spin 
model onto the continuum CFT. Moving along the H-N 
line we expect only a change in the non-universal quanti- 
ties like v, eoo and the values of the critical anisotropies 
themselves. We have seen that this is indeed the case, at 
least for these other two points: [A = 0.5, D c (0.5) = —1.2] 
(left part of fig. EJ) and [A = 1,D C (1) = -0.315] (both 
found by fixing A and varying D). For the former we 
get eoo = -2.00120, v = 2.44 and c = 0.5008 ± 0.0008, 
while for the latter we get v — 2.65, = —1.62651 and 
c = 0.498 ± 0.002. The corresponding gap exponents are 
estimated from the decay of /3l(— 1.2) and /3l(— 0.315), 
with the results v = 1.023 ± 0.009 and v = 1.003 ± 0.006, 
respectively. Notice that in fig. [futile two points have essen- 
tially the same /3-function. After all, eqs. I|35[) and (|37|l give 
(g c ) = C (1) L 1 /"a c /(2-v^ 1 ), where a c is the derivative 
of (g—g c ) with respect to the lattice parameter that is var- 
ied. Hence, if in fig.Elwe wish to compare the /3-functions 
of Ising type computed at fixed A with those computed 
at fixed D, we should multiply the latter by the (local) 
slope of H-N line d£>/dA ~ 1.697. We have checked that 
with a vertical shift of In 1.697 = 0.5289 the /3-function at 
(A = 1.185,1? = 0) collapses onto the other two, thereby 
indicating that moving along the line the c — 1/2 confor- 
mal structure is maintained. 

Concerning the discrepancy with ref. we should 
mention that a distinctive point of their analysis is the 
extrapolation for M — ► oo. On the other hand, we feel that 
the numerical procedure contains two weak points. First, 
the conclusion that v is not sufficiently close to 1, and 
the consequence that the effective dimensionality is not 
2, are drawn from the estimates of the correlation length 
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Fig. 8. Scaled gaps, divided by 2tt, plotted vs 1/L at the 
Ising transition (A = 1.1856, D — 0). Points represent the nu- 
merical values obtained with multi-target DMRG runs (L = 
20, 24, 28, 32, 36, 40, 44, 48, 56, 64, 72, 80 with M = 243) collect- 
ing eight excited states within M z = 0. Continuous lines the 
are best-fit whose intercepts are given in tabled together with 
the theoretical predictions of the scaling dimensions (the labels 
on the right indicate the multiplicities). 



Table 2. Spectrum of theoretical (dcFT = A + A+sec. indices) 
and numerical (see details in the text) conformal dimensions 
at the Ising transition of fig. |H| 



C?CFT 


CJnum 


Secondary 


Q 


[x multiplicity] 




Indices 




[xl] 




(0,0) 





1/8 [xl] 


0.1229 ± 0.0004 


(0,0) 


(or tt?) 


1 [xl] 


1 


(0,0) 


(or tt?) 


9/8 [x2] 


1.1185 ±0.0008 


(1,0) 


±2tt/L 




1.1218 ± 0.0008 


(0,1) 


=F2tt/L 


2 [x4] 


2.014 ±0.005 


(2,0) 


±4tt/L 




2.025 ±0.005 


(0,2) 


=F4tt/L 




2.018 ±0.005 


(1,0) 


±2tv/L 




2.018 ±0.005 


(0,1) 


=F2tt/L 



using a pure exponential law ignoring the contribution of 
algebraic prefactors. According to our experience, having 
a multi-target code at one's disposal, it would be better to 
read the values of v directly from the excitation gap or the 
scaling dimension. The second and more important point 
is the usage of the infinite-system algorithm. In sees. ^ 01 
and 0] we have given indications that this may be a risky 
procedure if one aims at very precise quantitative results, 
and that the finite-system algorithm should be preferred 
(see also refs. |2*lll9| regarding this question). 



C. Degli Esposti Boschi and F. Ortolani: Multi-target DMRG Methods for Quantum Phase Transitions 



11 



6 Concluding Remarks 

The aim of the this paper was to give a detailed expla- 
nation of the numerical method used to derive the results 
of ref. [17] and, more generally, to discuss how one can 
circumvent some of the problems of the DMRG close to 
criticality. The ideas are explicitly worked out taking some 
representative point in the A — D phase diagram for the 
Hamiltonian of eq. Q12JI. In sec. 01 we have argued that close 
to the transition lines from the Haldane phase (at least the 
ones with c = 1) the convergence is controlled by a char- 
acteristic length L* (M) that appears to be a consequence 
of the truncation to the M states with the largest weights 
in the block density matrix p\>. This fact is in line with 
similar results by other authors and this probably 

lies at the heart of the problem: in the TL the physical 
system behaves in a critical way but the DMRG proce- 
dure introduces a spurious length so that the numerical 
outcome is no longer scale invariant. 

The method that we are suggesting is based on the 
finite-system DMRG algorithm in order to reduce as much 
as possible the (environment 0) errors and has to be ap- 
plied in an intermediate range of L, not too small so that 
the signatures of scaling are visible but not too large as 
compared to L* (M). Then we make use of powerful finite- 
size scaling predictions, coming from CFT, to extract in- 
formation on the effective continuum model, like the "spin 
velocity" and the central charge, from the spectrum of low- 
lying excitations. The latter can be computed, within the 
DMRG, by taking pb as the average of the matrices asso- 
ciated with a given number of excited states that, in turn, 
are obtained with a Thick-restart variant of the Lanczos 
method |15| . 

We think that the methods discussed here should ap- 
ply also to the study of the critical behavior of other ID 
quantum systems, such as Hubbard models and their gen- 
eralizations |4*%>49I5()I51| . 

Finally, in sec. [S] we have reported some examples of 
c = 1/2 transitions between the Haldane and the Neel-like 
phase. In particular, we have re-examined and confirmed 
the 2D Ising nature of the critical point at D = |29) . 
which was the cause of a recent controversy 
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Appendix A: Thick-restart Lanczos Method 

Let us consider a n x n real symmetric matrix H repre- 
senting the Hamiltonian operator of a quantum system in 
a given orthonormal basis. We are interested in the low- 
est eigenvalues, {A}, and corresponding eigenvectors, x, 



defined by the equation ifx = Ax (A real, x £ M. n ). If 
the matrix H is large and only a small number of eigen- 
values are wanted, a projection-based method is generally 
used |f)2U53| . In these methods one usually builds orthogo- 
nal bases and then performs the Rayleigh-Ritz projection 
to extract the approximate solutions. When the matrix 
is symmetric, the Lanczos method is the most commonly 
used algorithm, taking advantage of the fact that the ac- 
tual matrix to diagonalize can be cast in tridiagonal form. 

Another method used for large matrices is the David- 
son one (and its variants) for which, however, the tridi- 
agonal form is not assured. When the Hamiltonian ma- 
trix is dominated by the diagonal elements, as occurs in 
quantum chemistry, the Davidson-Jacobi gives very good 
results This does not necessarily hold for a quantum 
spin system in the basis of eigenstates of local Sj and the 
preconditioner (especially during the initial steps of the 
infinite-system algorithm) can break some of the symme- 
tries of the Hamiltonian. Given that the choice of a partic- 
ular method depends on the physical system under study, 
in this appendix we focus on the characteristic points of 
the Lanczos method used for our numerical results. 

Given a starting vector ro <E W 1 , an orthogonal se- 
quence of vectors q; , i = 1,2, ... is generated recursively 
from the relations: 



qo = 
A) = ||ro|| 



0i-i 

Hqi - a;qi - /9i-iq l -i 
IN| , i = l,2,... 



(20) 



(21) 



The vectors q^, called Lanczos vectors, form an orthonor- 
mal basis for the Krylov subspaces: 



K m = span{r ,iiro, . . . ,H m 1 r } = span{qi, 
and satisfy the relations: 

Hqt = A-iq t -i + a»q 4 + Aqi+i , 



> } , 

(22) 



(23) 



which define a tridiagonal matrix, T m , as a representation 
of H in ICm : 



( Oil fa 

Pi cti 



V 



\ 



®"m— 1 Prn—1 
0m-l OL m / 



(24) 



The diagonalization of T m , with m — 1,2,..., gives eigen- 
values and eigenvectors, called Ritz values and Ritz vec- 
tors, as approximate solutions to the original problem. 
The computation of an eigenvalue with good accuracy 
often requires a large number of iterations and on most 
machines there is not enough memory to store all the 
Lanczos vectors. It is necessary to limit the number m 
of generated vectors and restart the procedure. Since the 
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algorithm defined by eqs. (121 )[) and H21(l can start with 
one vector ro, the usual way is to use the computed Ritz 
vector, if only one eigenvalue is wanted. If more than one 
eigenvalue is desired, we can freeze the converged ones and 
combine the other vectors into one starting element |53) . 
Other methods (Davidson-Liu or block Lanczos) are based 
on adding a block of states at each step but, of course, 
they require more computational resources. Typically, a 
restarting scheme needs significantly more iterations to 
compute a solution, but it saves memory usage. The algo- 
rithm developed by Wu and Simon |15| , summarized in the 
following, has the advantage of adding only one vector at 
every step, essentially preserving the tridiagonal structure 
of the projected matrix. 

The Thick-restart Lanczos method is based on the fol- 
lowing observation: let us consider an arbitrary element of 



(25) 



3=1 

from relation (|23[1 we have (setting ya = 0): 



Hq = (ft- 
i=i 

(26) 

so that the residual, r = y m f3 m qm+i, is always parallel to 
q m+ i independently from the vector q in IC m and has the 
same direction of the residual r m computed from eq. i|21|l 
and defining q m +i- This fact suggested to Wu and Simon 
a suitable restarting scheme. 

When restarting, one has first to determine an appro- 
priate number of Ritz vectors, say k (usually greater than 
the number of wanted eigenvalues): 



1. 



(27) 



corresponding to the Ritz values . The chosen Ritz vec- 
tors contain the best approximations available for the wan- 
ted eigenvectors. Since the matrix T m is symmetric, there 
is no reason to use any basis set different from its eigen- 
vectors: 

T m y (i) =\y (l) , i = l,...,k (28) 

(We denote with a tilde the quantities after the restart, to 
distinguish them from the corresponding ones before the 
restart.) Using these k vectors, the projected matrix of H, 
T k , is diagonal and the vectors satisfy the relations: 



Hqi = Xiqi + y$[3 m qk+i 



= 1, 



,,k 



(29) 



with cjfe + i = q m +i . So we can enlarge the basis {qi , . . . , qfc} 
with the vector q/c+i and, using the symmetry of the ma- 
trix H, the vector qfc+2 can be computed by the residual: 

k 

r fe+ i = Hq k+1 - & k+1 q k+1 - ^ 0mV$qj ■ (30) 

3=1 

Correspondingly, the matrix T k is extended by one row 
and one column into the matrix T k +\. The latter is not 



tridiagonal as the original T m but further steps follow the 
three-terms relations (1231) , defining on and (3i for i = k + 
2, fc+3, . . . , m, and at step m > k+ 1 we have the following 
form for T m : 



/Ai 



T — 



0i 

Pk °}k+l Pk+l 
Pk.+l Oik+2 



\ 



(31) 



with /3i = (3 m y 



(i) 



1, . . . , k. Since the vectors in the 



new subspace generated by {qi, . . . ,q m } have the same 
property that the residual is always parallel to q m+ i, the 
procedure can be repeatedly restarted, until we find a good 
convergence of the wanted eigenvalues. The matrix T m is 
no longer tridiagonal after the first restart, but it can still 
be stored and diagonalized in an efficient way. It is easy to 
arrange the algorithm so that q; and q\; occupy the same 
memory locations. 

We can summarize saying that the restarting scheme 
with k Ritz vectors q^ (we now drop the tilded notation), 
and a residual r k satisfying the rule: 

Hqi = Xiqi + ftqfe+i , (32) 

with qfc+i = rfc/||rfc||, is composed by an initialization: 

a k+i = {Qk+i,Hq k+1 ) 

r k+1 = Hq k+ i - a k+ iq k+1 - Pj<\ 3 (33) 

Pk+i = |kfe+i|| 

followed by typical Lanczos iterations (eq. Ij21|l ) for i = 
k + 2, k + 3, . . . , rn. 

This concludes our description of the Thick-restart 
Lanczos algorithm. We refer to the original paper for 
a detailed discussion of the errors and precision of the 
method. 



Appendix B: Brief Survey of FSS Theory 

Ideally, we should locate the critical points of a quan- 
tum system by looking at the smallest energy gap, AE, 
as a function of a (relevant) parameter g and identify g c 
through the condition AE(g c ) = 0. For an algebraic tran- 
sition the critical index v controls the opening of the gap, 
AE oc \g — g c \ u ■ However, when the critical point is ap- 
proached numerically we encounter two related problems. 
First, the system's size is necessarily finite and a true 
phase transition is forbidden. Second, if we are able to 
deal with sufficiently large systems close to g c , the energy 
gap may become so small as to be comparable with the 
errors introduced by the algorithm (or, ultimately, by the 
machine). Hence we need a prescription to infer the loca- 
tion of the critical point from nonzero values of AEl at 
finite L. 
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In FSS theory |2"4"ll45ll54| one usually invokes the fol- 
lowing ansatz (quantum Hamiltonian notations): 

AE L = jF{Q t C^L^\g~g c \, (34) 

so that the scaling variable C covers these two regimes 

(a) £ — > for g — > g c at fixed L. In this regime we should 
see a finite system at the infinite-size critical point. 
Since AEl is in fact an inverse correlation length £l, 
we expect the latter to be as large as possible, that 
is to say AEr, oc L~ l . Compatibility with eq. I|34|l 
then requires F(0) ^ 0. This is exactly the case for 
which the continuum description of CFT applies (sec. 

Hence, recalling that the scaling dimension of the 
gap-generating operator is (2 — v~ ), we may write 
F(0) = 2irv{2 - v~ Y ) and expand F(Q in a McLaurin 
series for ( -C 1: 

AE L = ^[(2-„- 1 ) + C^C + C^C 2 + ...] . (35) 
±j 

(b) £ — > oo, which means (L/^l) 1 ^ 1. This regime 
mimics the TL, in the sense that £l is large but finite 
because the system is slightly off-critical and L is suffi- 
ciently large so that scaling laws appear. Hence L must 
effectively cancel in eq. I)34|l . leaving just \g — g c \ v ■ This 
is possible if: 

F(Q ~ C C » 1 • (36) 

From (a) we derive the way to locate the critical points 
through FSS: Plot In (AEx) vs In L and look for the value 
of g that best gives a straight line with slope —1. Actually, 
we can be even more severe by looking for slope in the 
curves of the scaled gaps LAE^. As far as the index v 
is concerned, one may consider the finite-size /3-function, 
/3^ 1 (<7) = 9 In (AE^/dg, evaluated at g = g c as deter- 
mined above (apart from the sign): 

h(9c) = [F(C)/F'(C)] c =oi- 1/,y , (37) 

where F' denotes the derivative with respect to £ and 
F'(0) is also assumed to be nonzero. In principle /3l(<? c ) 
should vanish with an exponent \jv that represents the 
asymptotic slope in log-log scales. Alternatively, since the 
scaling region may be reachable only for very large L, one 
can calculate the discrete logarithmic derivative through 
a size increment L — > L + 5L: 

J_ ln/3 L+aL (ff c ) - In /3 L (g c ) 

v L ~ ln(L + 8L)-lnL ' 1 ' 

and this should converge to 1/v when L oo. 

Finally, the ansatz l|34|) can be generalized to other 
physical quantities that behave as Q(g) oc \g — g c \ ua near 
the critical point: 

Ql(9)=L-^F q (0 , (39) 

with the same scaling variable £ and scaling exponent zq. 
As above, in the critical regime (a) £ — > we shall require 



Fq(0) ^ 0, while in the off-critical regime (b) the scaling 
exponent vq emerges provided that Fq ~ Q" a for £ — > oo, 
and zq — vq/v. In many cases, Q 2 is a squared order 
parameter given by the asymptotic value of a certain cor- 
relation function, (Oq(O)Og(r)), that slightly away from 
the critical point behaves as: 

{O Q (0)O Q (r)) oc , (40) 

cIq being the scaling dimension of Oq. At this stage we 
can make the scaling argument of Ginsparg (sec. 5.1 of 
|lfi| ) evaluating eq. H4Q(I just at the correlation length it- 
self, r = £ cx \g - g c \- v thereby having (O a (0)O a (0> ~ 
Q 2 ~ Gg(l)|g — g c \ 2vda ■ Hence we find the scaling law 
cIq — vq/v, that tells us that zq in eq. (|39|l is nothing 
but the scaling dimension of the operator associated with 
the order parameter Q. In ref. ^7] we have used this prop- 
erty to derive the decay exponents of ordinary (transverse 
channel) and string (z-channel) correlation functions in 
c = 1 phases by applying eq. I|39() to the corresponding 
Neel and string order parameters evaluated at half chain. 
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